#ifndef COMB_UTIL_HPP__
#define COMB_UTIL_HPP__

#include <string>
#include <map>
#include <fstream>
#include <sstream>
#include <iostream>
#include <vector>
#include <cstdlib>
#include <cstring>
#include <valarray>

#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>

#include <boost/serialization/base_object.hpp>
#include <boost/serialization/vector.hpp>
#include <boost/serialization/map.hpp>
#include <boost/serialization/valarray.hpp>

#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/gzip.hpp>

#include <boost/multi_array.hpp>

#include "comb_types.hpp"


void error(const std::string &msg);

/** Give the reverse complement of a base.
 */
char rc(char c);

/** Return the reverse complement of a sequence.
 */
std::string rc(std::string seq);

/** Calculate the size of every sequence in a fasta file.
 */
std::map<std::string, unsigned int> build_fasta_index(const std::string &fname);

/** Open a PerM mapping file and figure out the type of problem, sequence/color.
 */
mapping_type get_mapping_type(const char *fname);

/** Open a PerM Mapping file and count the read length
 */
int get_read_length(const char *fname);

/** A class describing a position.
 */
class MapLoci
{
    public:
        MapLoci(std::string _id, unsigned int _pos, std::string _seq, char _strand);
        
        std::string     id;     // ID for current spot on reference
        unsigned int    pos;    // Position on genome where the match was made
        std::string     seq;    // Sequenced read from the reference genome
        char            strand; // Strand direction for this position
};

/** A read and the places it maps.
 */
class ColorRead 
{
    public:
        ColorRead(std::string _id, std::string _seq, std::vector<int> _qual, std::vector<MapLoci> _map_locs, int num);
        ColorRead();
        
        std::string id;                 // ID for current read
        std::string seq;                // Sequenced read
        std::vector<int> qual;          // Quality score for current read
        std::vector<MapLoci> map_locs;  // Contains MappingObj objects with data from reference genome
        unsigned int num;               // Number of reference position matches for this read
};


/** A read and the places it maps.
 */
class SeqRead 
{
    public:
        SeqRead(std::string _id, std::string _seq, std::vector<int> _qual, std::vector<MapLoci> _map_locs, int num);
        SeqRead();
        
        std::string id;                 // ID for current read
        std::string seq;                // Sequenced read
        std::vector<int> qual;          // Quality score for current read
        std::vector<MapLoci> map_locs;  // Contains MappingObj objects with data from reference genome
        unsigned int num;               // Number of reference position matches for this read
};

/** The record of a SNP candidate.
 */
class SnpRecord
{
    public:
        SnpRecord(unsigned int _id, std::string &_chr, unsigned int _start, unsigned int _end, std::valarray<float> &_subs, std::valarray<float> &_totl);
        SnpRecord();
        unsigned int    id;
        std::string     chr;
        unsigned int    start;
        unsigned int    end;
        std::string     seq;
        std::valarray<float> subs;
        std::valarray<float> totl;
};

/** What can be placed into a SNP candidate site.  
 */
class SnpModel
{
    public:
        SnpModel(double p[4], char c);
        SnpModel();
        double  prob[4];
        char    base;
};

/** Information used to call the SNP candidate.
 */
class ColorSnpCall
{
    public:
        ColorSnpCall(unsigned int _id, std::string &_chr, unsigned int _start, unsigned int _end, std::string _seq, std::valarray<float> &_subs, std::valarray<float> &_totl);
        ColorSnpCall();
        unsigned int    id;
        std::string     chr;
        unsigned int    start;
        unsigned int    end;
        std::string     seq;
        std::valarray<float> subs;
        std::valarray<float> totl;
        std::vector<std::vector<SnpModel> > blocks;
        std::vector<double> probs;
        int total_reads;
};


/** Information used to call the SNP candidate.
 */
class SeqSnpCall
{
    public:
        SeqSnpCall(unsigned int _id, std::string &_chr, unsigned int _start, unsigned int _end, std::string _seq, std::valarray<float> &_subs, std::valarray<float> &totl);
        SeqSnpCall();
        unsigned int    id;
        std::string     chr;
        unsigned int    start;
        unsigned int    end;
        std::string     seq;
        std::valarray<float> subs;
        std::valarray<float> totl;
        std::vector<SnpModel> blocks;
        std::vector<double> probs;
        int total_reads;
};

/** Split a string into a vector of strings using delim.
 */
void strtok(const std::string& str, const std::string &delim, std::vector<std::string>& tokens);

/** Read the fasta file.
 */
void parse_fasta(const char *fname, std::map<std::string, std::string> &seqs);

/** Write a fasta file to a stream.
 */
void write_fasta(std::ostream &os, std::map<std::string, std::string> &seqs);

/** Read the lengths of each sequence in a fasta file. 
 */
std::map<std::string, unsigned int> read_fasta_index(std::string fasta_index_fname);

// Global Variables.
extern std::map<char, std::map<char, char> > next_color;
extern std::map<char,int> iswap;
extern std::map<char, int> cswap;
extern std::map<char, std::map<char, int> > cswap2;
extern std::vector<SnpModel> haploid_models;
extern std::vector<std::vector<SnpModel> > diploid_models;

extern double haploid_priors_none[4][4];
extern double haploid_priors_high[4][4];
extern double haploid_priors_moderate[4][4];

extern const unsigned int SNP_CAND_FILE_TYPE;
extern const unsigned int CALL_RATES_FILE_TYPE; 
extern const unsigned int SUBS_TOTL_FILE_TYPE;
extern const unsigned int SNP_CALL_COL_FILE_TYPE;
extern const unsigned int SNP_CALL_SEQ_FILE_TYPE;
extern const unsigned int EXPRESSION_FILE_TYPE;

extern const unsigned int SNP_CAND_FILE_VERSION;
extern const unsigned int SNP_CALL_COL_FILE_VERSION;
extern const unsigned int SNP_CALL_SEQ_FILE_VERSION;
extern const unsigned int CALL_RATES_COL_FILE_VERSION;
extern const unsigned int CALL_RATES_SEQ_FILE_VERSION;
extern const unsigned int SUBS_TOTL_FILE_VERSION;
extern const unsigned int EXPRESSION_FILE_VERSION;

void init_globals();

std::string color_to_seq(std::string &seq);

typedef struct options_
{
    std::string map_file;
    std::string fasta_file;
    std::string cand_file;
    std::string call_file;
    std::string snp_file;
    std::string priors_file;
    std::string generic_file;
    std::string cons_file;
    std::vector<std::string> generic_file_list;

    bool unique;
    bool ignore_missing;

    bool only_show_changed; 

    bool has_cand;
    
    unsigned int max_map_loci;
    
    int consec_snps;

    std::vector<std::pair<int, float> > cover_threshold;
    
    unsigned int nbins;
    unsigned int read_count;

    bool diploid;
    unsigned int max_block_size;
    unsigned int old_bin_cnt;
    unsigned int new_bin_cnt;
    int normalize_freq;

    std::string priors_str;
    double haploid_priors[4][4];
    double diploid_priors[4][4];

    bool expr_all;
    bool expr_random;
    bool expr_unique;

    std::string map_count_file;

} prog_options;

/** Currently unused code for printing valarrays.  Useful in debugging.
 */
template<typename T>
void print_valarray(std::ostream& out, const std::valarray<T>& a)
{
    out << "{ ";
    for(size_t i = 0; i < a.size(); ++i)
        out << " " << a[i];
    out << " }";
}

/** Currently unused code for printing valarrays.  Useful in debugging.
 */
/*
template< template< typename T> class U, typename T>
std::ostream &operator<<(std::ostream &out, const U<T>& x)
{
    print_valarray(out, static_cast<std::valarray<T> >(x));
    return out;
}
*/
#endif
